home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / docs / apps / circuits / spice2g6.z / spice2g6 / spice / Fortran / jfet.f < prev    next >
Encoding:
Text File  |  1989-02-03  |  10.2 KB  |  354 lines
  1.       subroutine jfet
  2.       implicit double precision (a-h,o-z)
  3. c
  4. c     this routine processes jfets for dc and transient analyses.
  5. c
  6. c spice version 2g.6  sccsid=tabinf 3/15/83
  7.       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
  8.      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
  9.      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
  10.      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
  11.      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
  12.      5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval,
  13.      6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt,
  14.      7   irowno,jcolno,nttbr,nttar,lvntmp
  15. c spice version 2g.6  sccsid=cirdat 3/15/83
  16.       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
  17.      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc
  18. c spice version 2g.6  sccsid=status 3/15/83
  19.       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
  20.      1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
  21.      2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
  22. c spice version 2g.6  sccsid=knstnt 3/15/83
  23.       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
  24.      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
  25.      2   pivtol,pivrel
  26. c spice version 2g.6  sccsid=blank 3/15/83
  27.       common /blank/ value(200000)
  28.       integer nodplc(64)
  29.       complex cvalue(32)
  30.       equivalence (value(1),nodplc(1),cvalue(1))
  31. c
  32. c
  33.       dimension vgso(1),vgdo(1),cgo(1),cdo(1),cgdo(1),gmo(1),gdso(1),
  34.      1   ggso(1),ggdo(1),qgs(1),cqgs(1),qgd(1),cqgd(1)
  35.       equivalence (vgso(1),value( 1)),(vgdo(1),value( 2)),
  36.      1            (cgo (1),value( 3)),(cdo (1),value( 4)),
  37.      2            (cgdo(1),value( 5)),(gmo (1),value( 6)),
  38.      3            (gdso(1),value( 7)),(ggso(1),value( 8)),
  39.      4            (ggdo(1),value( 9)),(qgs (1),value(10)),
  40.      5            (cqgs(1),value(11)),(qgd (1),value(12)),
  41.      6            (cqgd(1),value(13))
  42. c
  43. c
  44.       loc=locate(13)
  45.    10 if ((loc.eq.0).or.(nodplc(loc+25).ne.0)) return
  46.       locv=nodplc(loc+1)
  47.       node1=nodplc(loc+2)
  48.       node2=nodplc(loc+3)
  49.       node3=nodplc(loc+4)
  50.       node4=nodplc(loc+5)
  51.       node5=nodplc(loc+6)
  52.       locm=nodplc(loc+7)
  53.       ioff=nodplc(loc+8)
  54.       type=nodplc(locm+2)
  55.       locm=nodplc(locm+1)
  56.       loct=nodplc(loc+19)
  57. c
  58. c  dc model parameters
  59. c
  60.       area=value(locv+1)
  61.       vto=value(locm+1)
  62.       beta=value(locm+2)*area
  63.       xlamb=value(locm+3)
  64.       gdpr=value(locm+4)*area
  65.       gspr=value(locm+5)*area
  66.       csat=value(locm+9)*area
  67.       vcrit=value(locm+16)
  68. c
  69. c  initialization
  70. c
  71.       icheck=1
  72.       go to (100,20,30,50,60,70), initf
  73.    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
  74.       vds=type*value(locv+2)
  75.       vgs=type*value(locv+3)
  76.       vgd=vgs-vds
  77.       go to 300
  78.    25 if(ioff.ne.0) go to 40
  79.       vgs=-1.0d0
  80.       vgd=-1.0d0
  81.       go to 300
  82.    30 if (ioff.eq.0) go to 100
  83.    40 vgs=0.0d0
  84.       vgd=0.0d0
  85.       go to 300
  86.    50 vgs=vgso(lx0+loct)
  87.       vgd=vgdo(lx0+loct)
  88.       go to 300
  89.    60 vgs=vgso(lx1+loct)
  90.       vgd=vgdo(lx1+loct)
  91.       go to 300
  92.    70 xfact=delta/delold(2)
  93.       vgso(lx0+loct)=vgso(lx1+loct)
  94.       vgs=(1.0d0+xfact)*vgso(lx1+loct)-xfact*vgso(lx2+loct)
  95.       vgdo(lx0+loct)=vgdo(lx1+loct)
  96.       vgd=(1.0d0+xfact)*vgdo(lx1+loct)-xfact*vgdo(lx2+loct)
  97.       cgo(lx0+loct)=cgo(lx1+loct)
  98.       cdo(lx0+loct)=cdo(lx1+loct)
  99.       cgdo(lx0+loct)=cgdo(lx1+loct)
  100.       gmo(lx0+loct)=gmo(lx1+loct)
  101.       gdso(lx0+loct)=gdso(lx1+loct)
  102.       ggso(lx0+loct)=ggso(lx1+loct)
  103.       ggdo(lx0+loct)=ggdo(lx1+loct)
  104.       go to 110
  105. c
  106. c  compute new nonlinear branch voltages
  107. c
  108.   100 vgs=type*(value(lvnim1+node2)-value(lvnim1+node5))
  109.       vgd=type*(value(lvnim1+node2)-value(lvnim1+node4))
  110.   110 delvgs=vgs-vgso(lx0+loct)
  111.       delvgd=vgd-vgdo(lx0+loct)
  112.       delvds=delvgs-delvgd
  113.       cghat=cgo(lx0+loct)+ggdo(lx0+loct)*delvgd+ggso(lx0+loct)*delvgs
  114.       cdhat=cdo(lx0+loct)+gmo(lx0+loct)*delvgs+gdso(lx0+loct)*delvds
  115.      1   -ggdo(lx0+loct)*delvgd
  116. c
  117. c  bypass if solution has not changed
  118. c
  119.       if (initf.eq.6) go to 200
  120.       tol=reltol*dmax1(dabs(vgs),dabs(vgso(lx0+loct)))+vntol
  121.       if (dabs(delvgs).ge.tol) go to 200
  122.       tol=reltol*dmax1(dabs(vgd),dabs(vgdo(lx0+loct)))+vntol
  123.       if (dabs(delvgd).ge.tol) go to 200
  124.       tol=reltol*dmax1(dabs(cghat),dabs(cgo(lx0+loct)))+abstol
  125.       if (dabs(cghat-cgo(lx0+loct)).ge.tol) go to 200
  126.       tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
  127.       if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
  128.       vgs=vgso(lx0+loct)
  129.       vgd=vgdo(lx0+loct)
  130.       vds=vgs-vgd
  131.       cg=cgo(lx0+loct)
  132.       cd=cdo(lx0+loct)
  133.       cgd=cgdo(lx0+loct)
  134.       gm=gmo(lx0+loct)
  135.       gds=gdso(lx0+loct)
  136.       ggs=ggso(lx0+loct)
  137.       ggd=ggdo(lx0+loct)
  138.       go to 900
  139. c
  140. c  limit nonlinear branch voltages
  141. c
  142.   200 ichk1=1
  143.       call pnjlim(vgs,vgso(lx0+loct),vt,vcrit,icheck)
  144.       call pnjlim(vgd,vgdo(lx0+loct),vt,vcrit,ichk1)
  145.       if (ichk1.eq.1) icheck=1
  146.       call fetlim(vgs,vgso(lx0+loct),vto)
  147.       call fetlim(vgd,vgdo(lx0+loct),vto)
  148. c
  149. c  determine dc current and derivatives
  150. c
  151.   300 vds=vgs-vgd
  152.       if (vgs.gt.-5.0d0*vt) go to 310
  153.       ggs=-csat/vgs+gmin
  154.       cg=ggs*vgs
  155.       go to 320
  156.   310 evgs=dexp(vgs/vt)
  157.       ggs=csat*evgs/vt+gmin
  158.       cg=csat*(evgs-1.0d0)+gmin*vgs
  159.   320 if (vgd.gt.-5.0d0*vt) go to 330
  160.       ggd=-csat/vgd+gmin
  161.       cgd=ggd*vgd
  162.       go to 340
  163.   330 evgd=dexp(vgd/vt)
  164.       ggd=csat*evgd/vt+gmin
  165.       cgd=csat*(evgd-1.0d0)+gmin*vgd
  166.   340 cg=cg+cgd
  167. c
  168. c  compute drain current and derivitives for normal mode
  169. c
  170.   400 if (vds.lt.0.0d0) go to 450
  171.       vgst=vgs-vto
  172. c
  173. c  normal mode, cutoff region
  174. c
  175.       if (vgst.gt.0.0d0) go to 410
  176.       cdrain=0.0d0
  177.       gm=0.0d0
  178.       gds=0.0d0
  179.       go to 490
  180. c
  181. c  normal mode, saturation region
  182. c
  183.   410 betap=beta*(1.0d0+xlamb*vds)
  184.       twob=betap+betap
  185.       if (vgst.gt.vds) go to 420
  186.       cdrain=betap*vgst*vgst
  187.       gm=twob*vgst
  188.       gds=xlamb*beta*vgst*vgst
  189.       go to 490
  190. c
  191. c  normal mode, linear region
  192. c
  193.   420 cdrain=betap*vds*(vgst+vgst-vds)
  194.       gm=twob*vds
  195.       gds=twob*(vgst-vds)+xlamb*beta*vds*(vgst+vgst-vds)
  196.       go to 490
  197. c
  198. c  compute drain current and derivitives for inverse mode
  199. c
  200.   450 vgdt=vgd-vto
  201. c
  202. c  inverse mode, cutoff region
  203. c
  204.       if (vgdt.gt.0.0d0) go to 460
  205.       cdrain=0.0d0
  206.       gm=0.0d0
  207.       gds=0.0d0
  208.       go to 490
  209. c
  210. c  inverse mode, saturation region
  211. c
  212.   460 betap=beta*(1.0d0-xlamb*vds)
  213.       twob=betap+betap
  214.       if (vgdt.gt.-vds) go to 470
  215.       cdrain=-betap*vgdt*vgdt
  216.       gm=-twob*vgdt
  217.       gds=xlamb*beta*vgdt*vgdt-gm
  218.       go to 490
  219. c
  220. c  inverse mode, linear region
  221. c
  222.   470 cdrain=betap*vds*(vgdt+vgdt+vds)
  223.       gm=twob*vds
  224.       gds=twob*vgdt-xlamb*beta*vds*(vgdt+vgdt+vds)
  225. c
  226. c  compute equivalent drain current source
  227. c
  228.   490 cd=cdrain-cgd
  229.       if (mode.ne.1) go to 500
  230.       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
  231.       if (initf.eq.4) go to 500
  232.       go to 700
  233. c
  234. c  charge storage elements
  235. c
  236.   500 czgs=value(locm+6)*area
  237.       czgd=value(locm+7)*area
  238.       phib=value(locm+8)
  239.       twop=phib+phib
  240.       fcpb=value(locm+12)
  241.       fcpb2=fcpb*fcpb
  242.       f1=value(locm+13)
  243.       f2=value(locm+14)
  244.       f3=value(locm+15)
  245.       czgsf2=czgs/f2
  246.       czgdf2=czgd/f2
  247.       if (vgs.ge.fcpb) go to 510
  248.       sarg=dsqrt(1.0d0-vgs/phib)
  249.       qgs(lx0+loct)=twop*czgs*(1.0d0-sarg)
  250.       capgs=czgs/sarg
  251.       go to 520
  252.   510 qgs(lx0+loct)=czgs*f1+czgsf2*(f3*(vgs-fcpb)
  253.      1   +(vgs*vgs-fcpb2)/(twop+twop))
  254.       capgs=czgsf2*(f3+vgs/twop)
  255.   520 if (vgd.ge.fcpb) go to 530
  256.       sarg=dsqrt(1.0d0-vgd/phib)
  257.       qgd(lx0+loct)=twop*czgd*(1.0d0-sarg)
  258.       capgd=czgd/sarg
  259.       go to 560
  260.   530 qgd(lx0+loct)=czgd*f1+czgdf2*(f3*(vgd-fcpb)
  261.      1   +(vgd*vgd-fcpb2)/(twop+twop))
  262.       capgd=czgdf2*(f3+vgd/twop)
  263. c
  264. c  store small-signal parameters
  265. c
  266.   560 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
  267.       if (initf.ne.4) go to 600
  268.       value(lx0+loct+9)=capgs
  269.       value(lx0+loct+11)=capgd
  270.       go to 1000
  271. c
  272. c  transient analysis
  273. c
  274.   600 if (initf.ne.5) go to 610
  275.       qgs(lx1+loct)=qgs(lx0+loct)
  276.       qgd(lx1+loct)=qgd(lx0+loct)
  277.   610 call intgr8(geq,ceq,capgs,loct+9)
  278.       ggs=ggs+geq
  279.       cg=cg+cqgs(lx0+loct)
  280.       call intgr8(geq,ceq,capgd,loct+11)
  281.       ggd=ggd+geq
  282.       cg=cg+cqgd(lx0+loct)
  283.       cd=cd-cqgd(lx0+loct)
  284.       cgd=cgd+cqgd(lx0+loct)
  285.       if (initf.ne.5) go to 700
  286.       cqgs(lx1+loct)=cqgs(lx0+loct)
  287.       cqgd(lx1+loct)=cqgd(lx0+loct)
  288. c
  289. c  check convergence
  290. c
  291.   700 if (initf.ne.3) go to 710
  292.       if (ioff.eq.0) go to 710
  293.       go to 750
  294.   710 if (icheck.eq.1) go to 720
  295.       tol=reltol*dmax1(dabs(cghat),dabs(cg))+abstol
  296.       if (dabs(cghat-cg).ge.tol) go to 720
  297.       tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol
  298.       if (dabs(cdhat-cd).le.tol) go to 750
  299.   720 noncon=noncon+1
  300.   750 vgso(lx0+loct)=vgs
  301.       vgdo(lx0+loct)=vgd
  302.       cgo(lx0+loct)=cg
  303.       cdo(lx0+loct)=cd
  304.       cgdo(lx0+loct)=cgd
  305.       gmo(lx0+loct)=gm
  306.       gdso(lx0+loct)=gds
  307.       ggso(lx0+loct)=ggs
  308.       ggdo(lx0+loct)=ggd
  309. c
  310. c  load current vector
  311. c
  312.   900 ceqgd=type*(cgd-ggd*vgd)
  313.       ceqgs=type*((cg-cgd)-ggs*vgs)
  314.       cdreq=type*((cd+cgd)-gds*vds-gm*vgs)
  315.       value(lvn+node2)=value(lvn+node2)-ceqgs-ceqgd
  316.       value(lvn+node4)=value(lvn+node4)-cdreq+ceqgd
  317.       value(lvn+node5)=value(lvn+node5)+cdreq+ceqgs
  318. c
  319. c  load y matrix
  320. c
  321.       locy=lvn+nodplc(loc+20)
  322.       value(locy)=value(locy)+gdpr
  323.       locy=lvn+nodplc(loc+21)
  324.       value(locy)=value(locy)+ggd+ggs
  325.       locy=lvn+nodplc(loc+22)
  326.       value(locy)=value(locy)+gspr
  327.       locy=lvn+nodplc(loc+23)
  328.       value(locy)=value(locy)+gdpr+gds+ggd
  329.       locy=lvn+nodplc(loc+24)
  330.       value(locy)=value(locy)+gspr+gds+gm+ggs
  331.       locy=lvn+nodplc(loc+9)
  332.       value(locy)=value(locy)-gdpr
  333.       locy=lvn+nodplc(loc+10)
  334.       value(locy)=value(locy)-ggd
  335.       locy=lvn+nodplc(loc+11)
  336.       value(locy)=value(locy)-ggs
  337.       locy=lvn+nodplc(loc+12)
  338.       value(locy)=value(locy)-gspr
  339.       locy=lvn+nodplc(loc+13)
  340.       value(locy)=value(locy)-gdpr
  341.       locy=lvn+nodplc(loc+14)
  342.       value(locy)=value(locy)+gm-ggd
  343.       locy=lvn+nodplc(loc+15)
  344.       value(locy)=value(locy)-gds-gm
  345.       locy=lvn+nodplc(loc+16)
  346.       value(locy)=value(locy)-ggs-gm
  347.       locy=lvn+nodplc(loc+17)
  348.       value(locy)=value(locy)-gspr
  349.       locy=lvn+nodplc(loc+18)
  350.       value(locy)=value(locy)-gds
  351.  1000 loc=nodplc(loc)
  352.       go to 10
  353.       end
  354.